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In this work, the solvation and electronic structure of the aqueous chloride ion solution was 
investigated using Density Functional Theory (DFT) based ab initio molecular dynamics (AIMD). 

From an analysis of radial distribution functions, coordination numbers, and solvation structures, 
we found that exact exchange (Axx) and non-local van der Waals (vdW) interactions effectively 
weaken the interactions between the Cl~ ion and the first solvation shell. With a Cl-O coordination 
number in excellent agreement with experiment, we found that most configurations generated with 
vdW-inclusive hybrid DFT exhibit 6-fold coordinated distorted trigonal prism structures, which is 
indicative of a significantly disordered first solvation shell. By performing a series of band structure 
calculations on configurations generated from AIMD simulations with varying DFT potentials, we 
found that the solvated ion orbital energy levels (unlike the band structure of liquid water) strongly 
depend on the underlying molecular structures. In addition, these orbital energy levels were also 
significantly affected by the DFT functional employed for the electronic structure; as the fraction of 
Axx was increased, the gap between the highest occupied molecular orbital of Cl“ and the valence 
band maximum of liquid water steadily increased towards the experimental value. 

PACS numbers: 


I. INTRODUCTION 

The nature of the interaction between the hydrogen 
bond (HB) network of liquid water and the solvated chlo¬ 
ride ion (a member of the Hofmeister series) is currently 
a topic under intense research due to its fundamental 
importance in biochemistry, atmospheric chemistry, and 
geological processes.Experimentally, the interaction 
between the solvated chloride ion, Cl“, and the surround¬ 
ing aqueous environment has been successfully probed 
by a variety of techniques.In particular, X-ray and 
neutron scattering as well as X-ray absorption measure¬ 
ments, have been used to investigate the solvation struc¬ 
ture of aqueous chloride solutions, providing experimen¬ 
tal determination of the Cl-O and Cl-H radial distribu¬ 
tion functions (RDF), gci-oif) and 5 ci-h(?'), and the as¬ 
sociated coordination number, i.e., the number of water 
molecules populating the first solvation shell surrounding 
a given Cl“ ion.®“^® In addition, the electronic structure 
associated with this fundamental aqueous ionic solution 
was studied by Delahay^^ and Winter et in which 

state-of-the-art photoemission spectroscopy (PES) was 
utilized to characterize the valence 3p bands of the sol¬ 
vated Cl“ ion, which were found to be approximately 
1.25-1.50 eV above the valence band maximum (VBM) 
of liquid water. 

From a theoretical and computational point of view, 
first-principles based computational methods such as ab 
initio molecular dynamics (AIMD) have become power¬ 
ful tools in the study of condensed-phase systems such 
as liquids^®’^^ and aqueous ionic solutions.With the 
AIMD technique, the nuclear potential energy surface is 


generated “on the fly” from the electronic ground state^° 
without the need for empirical input, thereby allowing for 
a quantum mechanical treatment of not only the struc¬ 
ture and dynamics of a given molecular system of inter¬ 
est, but also its electronic and dielectric properties, as 
well as potential chemical reactions^ {i.e., the break¬ 
ing and forming of chemical bonds). Since the initial 
pioneering simulations of liquid water,AIMD has 
been applied to many complex problems in biology, chem¬ 
istry, and energy research, e.g., the designing of efficient 
catalysts for hydrogen production^^ and the modeling of 
the auto-ionization^^ and electrocatalytic splitting of wa¬ 
ter,to name a few. 

In particular, AIMD simulations employing density 
functional theory (DFT) as the source of the underly¬ 
ing quantum mechanical potential have been used to 
generate reasonably accurate microscopic descriptions of 
the structure of the aqueous chloride ion solution, 
with RDFs and associated coordination numbers that 
were consistent with the available experimental data. 
These studies found that the predicted solvation struc¬ 
ture in the aqueous chloride ion solution is moderately 
dependent on the exchange-correlation (XC) potential, 

with semi-local generalized-gradient approximation^^ds 

(GGA) functionals such as PBE^^ BLYP^^ds yield¬ 
ing similar predictions to hybrid functionals such as 
PBEO.^® In addition, the HB network of water was found 
to be only locally perturbed by the presence of the sol¬ 
vated ion, a finding which is consistent with the general 
expectation that Gl“ is a weak disrupter of the aqueous 
environment.However, these studies also demon¬ 
strated that the electronic structure of this aqueous ionic 
solution is very sensitive to the adopted XC approxima- 
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In particular, higher-level electronic structure 
calculations performed on GGA-DFT generated config¬ 
urations (structures) were found to significantly under¬ 
estimate the energy levels of the solvated Gl“ ion, with 
resulting energetics that can even be qualitatively incor¬ 
rect.®'^’®® 

In the study of an aqueous ionic solution, an accurate 
theoretical description of the underlying HB network is 
a prerequisite for characterizing and understanding the 
interaction between a solvated or embedded ion and its 
surrounding environment. In neat liquid water, the pre¬ 
dictive power of DFT-based AIMD in the microscopic 
description of the HB network critically depends on the 
accuracy of the underlying XC functional utilized in the 
quantum mechanical treatment of the electronic degrees 
of freedom. In this regard, it is now clear that the most 
widely used XC functionals, i.e., the class of functionals 
based on the GGA, have severe limitations when applied 
to liquid water®®^^® as well as the crystalline phases of 
ice.®®“®^ Most notably, GGA-DFT suffers from the pres¬ 
ence of self-interaction error (SIE)®®’®® and the neglect of 
non-local electron correlation effects that are responsible 
for van der Waals (vdW) or dispersion interactions; as a 
result, these deficiencies in the GGA XC potential man¬ 
ifest as significant overstructuring and excessively slug¬ 
gish dynamics in aqueous systems such as ambient liq¬ 
uid water. Beyond the choice of the XC functional, most 
AIMD simulations of liquid water and aqueous ionic solu¬ 
tions performed to date have adopted classical mechanics 
for the nuclear equations of motion and therefore com¬ 
pletely neglect nuclear quantum effects (NQE)—another 
approximation that has been deemed insufficient for a 
highly accurate quantitative description of the micro¬ 
scopic structure and HB network in aqueous systems. In 
the case of liquid water, light atoms such as hydrogen de¬ 
viate significantly from classical behavior even at room 
temperature,as evidenced by experimental isotope 
effect studies which demonstrated a softening of the liq¬ 
uid structure (i.e., in the comparison of H 2 O to D 2 O).®® 
Hence the neglect of NQE in aqueous systems such as 
ambient liquid water leads to overstructuring in the pre¬ 
dicted RDFs.®^’®''’^®® 

A commonly adopted method to alleviate the delete¬ 
rious effects of SIE in GGA-DFT is the use of hybrid 
XG functionals, wherein a fraction of exact (or Hartree- 
Fock) exchange {E^x) is included in the density func¬ 
tional approximation. Due to the relatively high compu¬ 
tational cost associated with these XC functionals, ap¬ 
plications of hybrid DFT have mostly been restricted 
to small gas-phase clusters of water,although re¬ 
cently hybrid functionals have been applied in the study 
of several crystalline phases of ice®^“®"^’^®^ and liquid wa- 
ter. 50,51,64, 7 ®,79, 10 ® comparison to GGAs, these studies 
demonstrated that the energetic, structural, and vibra¬ 
tional properties of these aqueous systems, as predicted 
by hybrid DFT calculations, are generally in closer agree¬ 
ment with the available experimental data.®^“®^’®®“^®® In¬ 
deed, an accurate microscopic description of the HB net¬ 


work by hybrid XC functionals has a non-negligible ef¬ 
fect on the theoretical characterization of the solvated 
chloride ion. In this regard, it was found that the use 
of hybrid DFT was crucial in obtaining a qualitatively 
correct energy difference between Cl“ and the VBM of 
liquid water. In addition, both hybrid and GGA XC 
functionals lack the ability to describe vdW/dispersion 
interactions, which arise from non-local dynamical elec¬ 
tron correlation and have a substantial effect on the mi¬ 
croscopic structure of condensed-phase aqueous systems. 
In fact, the explicit inclusion of vdW interactions in DFT 
has been shown to significantly improve upon the theo¬ 
retical description of the transition pressures among the 
high-pressure phases of ice®^’®'*’’®® and the predicted equi¬ 
librium density of liquid water. ®^’^^’^® While many recent 
studies have concluded that the structure of liquid wa¬ 
ter significantly softens when vdW interactions are ac¬ 
counted for in the underlying XC potential, the extent to 
which these non-local forces affect the structure of liquid 
water is largely dependent upon the given approach uti¬ 
lized to facilitate vdW-inclusive DFT.®^’™“^®’^®® Hence, 
the inclusion of vdW or dispersion interactions in the un¬ 
derlying XC potential is expected to have a non-negligible 
effect on the theoretical description of the HB network 
and must be accounted for in the study of aqueous ionic 
solutions such as the solvated chloride ion. 

Recently it was shown that utilization of the hy¬ 
brid PBEO functional,which includes 25% ex¬ 
act exchange, in conjunction with a fully self- 
consistent (SC) implementation of the density-dependent 
vdW/dispersion correction of Tkatchenko and Schef- 
flerio® (TS-vdW), i.e., the PBEO-bTS-vdW(SC) XC func¬ 
tional, yields an oxygen-oxygen structure factor, Soo{Q), 
and corresponding RDF, 500 (^)i that are in quanti¬ 
tative agreement with the best available experimental 
data.^®®’^®® This level of agreement between ab initio 
simulations and experiment was attributed to an increase 
in the relative population of water molecules in the in¬ 
terstitial region {i.e., the region between the first and 
second coordination shells), a collective reorganization 
in the liquid phase which is facilitated by a weakening of 
the HB strength by the use of a hybrid XC functional, 
coupled with a relative stabilization of the resultant dis¬ 
ordered liquid water configurations by the inclusion of 
long-range vdW/dispersion interactions. In fact, this in¬ 
creasingly more accurate description of the underlying 
HB network in liquid water also yielded other correlation 
functions, such as the oxygen-hydrogen RDF, 5 oh( 7'), 
and the higher-order oxygen-oxygen-oxygen triplet an¬ 
gular distribution, Pooo{(^), which encodes the degree 
of local tetrahedrality, as well as electrostatic properties, 
such as the effective molecular dipole moment, that are 
in much better agreement with experiment.^®® In this 
regard, the overall agreement between experiment and 
the PBEO-l-TS-vdW(SC) description of the microscopic 
structure of ambient liquid water is indeed a very promis¬ 
ing starting point for accurately investigating the struc¬ 
tural and energetic properties that characterize the aque- 
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ous chloride ion solution, which is the main focus of the 
work reported herein. 


II. COMPUTATIONAL DETAILS 

In this work, we have systematically performed a se¬ 
ries of Car-Parrinello AIMD simulations^*^ of the aqueous 
chloride ion solution at ambient conditions using a hier¬ 
archy of different XC functionals. The sequence of XC 
functionals employed herein includes the standard semi¬ 
local GGA of Perdew, Burke, and Ernzerhof (PBE),^® 
the corresponding hybrid PBEO^®’*^**^ which includes 25% 
exact exchange, and the self-consistent (SG) dispersion- 
corrected analogs^**^ thereof, i.e., PBE-l-TS-vdW(SC) 
and PBEO-l-TS-vdW(SG), based on the Tkatchenko- 
Scheffler*^*^® density-dependent vdW/dispersion func¬ 
tional. 

All of these AIMD simulations were performed in the 
canonical (NVT) ensemble using periodic simple cubic 
simulation cells containing Cl“ ion in (H20)63 with lat¬ 
tice parameters set to reproduce the experimental density 
of liquid water at ambient conditions. All of the AIMD 
simulations were initially equilibrated for approximately 
6 ps and then continued for an additional 20-50 ps for 
data collection (with the AIMD simulations employing 
the hybrid PBEO functional at the lower end of this range 
of simulation times). Since a classical treatment of the 
nuclear degrees of freedom was found to be insufficient for 
a quantitatively accurate description of the microscopic 
structure of ambient liquid water, we have performed all 
AIMD simulations in this work at the elevated tempera¬ 
ture of 330 K, a technique that is suggested by the lowest- 
order perturbative expansion of the free energy (in h) to 
account for the quantum mechanical nature of the nu¬ 
clear degrees of freedom. ^*^® In practice, this increase of 
approximately 30 K in the simulation temperature has 
been found to mimic the nuclear quantum effects (NQE) 
in structural quantities such as the oxygen-oxygen radial 
distribution function (500 (''’)) hi both DET®^ and force 
field®®’®® based MD simulations of liquid water. 

All calculations reported herein were performed within 
the plane-wave and pseudopotential framework and uti¬ 
lized a modified development version of the Quantum 
ESPRESSO (QE) software package.^*®® To meet the ad¬ 
ditional computational demands associated with large- 
scale AIMD simulations based on hybrid XG function¬ 
als, we have employed a linear scaling 0{N) exact ex¬ 
change algorithm that exploits the natural sparsity as¬ 
sociated with the real-space maximally localized Wan- 
nier function (MLWE)^*^® representation of the occupied 
Kohn-Sham electronic states.In addition, we have 
also developed and utilized a linear scaling 0{N) self- 
consistent implementation of the TS-vdW dispersion cor¬ 
rection,*^®^ which provides a framework for computing 
atomic Cg dispersion coefficients as explicit functionals 
of the charge density, i.e., Cq^ab = Cq^ab[p{'^)]i thereby 
accounting for the local chemical environment surround¬ 


ing each atom.*®® More explicit descriptions of the the¬ 
oretical methods employed herein can be found in Ref. 
[103]. 

The Gar-Parrinello (GP)^® equations of motion for the 
nuclear and electronic degrees of freedom were integrated 
using the standard Verlet algorithm and a time step of 
4.0 a.u. (ss 0.1 fs). To ensure an adiabatic separa¬ 
tion between the electronic and nuclear degrees of free¬ 
dom in the CP dynamics, we used a fictitious electronic 
mass of 300 a.u., which was found to be a reasonable 
choice for the simulation of water,®® and the nuclear 
mass of deuterium for each hydrogen atom. All elec¬ 
tronic wavefunctions were expanded using a plane wave 
basis set with a kinetic energy cutoff of 72 Ry, and the 
interactions between the valence electrons and the ions 
(consisting of the nuclei and their corresponding frozen- 
core electrons) were treated with Troullier-Martins type 
norm-conserving pseudopotentials.**® Ionic temperatures 
were controlled with massive Nose-Hoover chain ther¬ 
mostats,**^’**® each with a chain length of 4.**® A neu¬ 
tralizing background charge was included in the Ewald 
summation of the electrostatic energy to compensate for 
the net negative charge of the chloride ion solution. 

To examine the electronic structure surrounding the 
solvated chloride ion, we have performed wavefunction 
optimizations (i.e., self-consistent solutions of the non¬ 
linear Kohn-Sham equations) using the PBE, PBEO (25% 
exact exchange), and BHLYP**'’ (50% exact exchange) 
XG functionals on a set of configurations (structures) 
selected at even intervals along a given MD trajectory. 
Real-space integration of the orbital densities surround¬ 
ing the solvated chloride ion was performed within a 
sphere of radius i? = 1.5 A to determine the ionic contri¬ 
bution to the overall electronic structure and to identify 
the highest-occupied molecular orbital (HOMO) states. 
All PBE and PBEO wavefunction optimizations were per¬ 
formed with the QE software package.*®® Galculations 
with the BHLYP XG functional were performed with the 
CP2K suite of programs**® using a split-valence triple-^ 
basis set (appended with an additional set of polarization 
functions) in conjunction with Goedecker-Teter-Hutter 
type pseudopotentials**® and a density cut-off of 400Ry. 

III. RESULTS AND DISCUSSIONS 

A. Solvation Structure of the Aqueous Cl“ Ion 

We begin our investigation of the aqueous chloride 
ion solution by considering the chloride-oxygen and 
chloride-hydrogen RDEs, gc\-oii^) and ffci-H(''’); as 
shown in Fig. 1. The numerical values for the intensi¬ 
ties and h™j“jj) and positions and 

of the first maximum, the positions (r^^Q and r™“jj) of 
the first minimum, as well as the corresponding coordina¬ 
tion numbers (nci-o and nci-n) computed from these 
RDFs are provided in Table I along with the available 
experimental findings. Here, the nci-o and nqi-H asso- 
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TABLE I: Numerical values for the intensities (hcn?o s-nd hcfTn) a-nd positions and rcilTH ia of tho first maximum 

in the RDFs depicted in Fig. 1, the positions (rci‘-o and cciIIh in A) of the first minimum in the RDFs depicted in Fig. 1, as 
well as the corresponding coordination numbers (nci-o and nci-n). Available experimental findings for these quantities have 
been compiled in the bottom rows. 


Method 

I, max _max ^min „ 

^Cl-O T-ci-O T-ci-O Tlci-O 

Lmax _max „ 

«C1-H 7-ci h 7-ci h TIcI-H 

PBE 

PBE-bTS-vdW(SC) 

PBEO 

PBE0“’*' 

PBEO-tTS-vdW(SC) 

Expt.'^ 

Exp!.'* 

2.67T0.11 3.11T0.02 3.64T0.03 5.5±0.2 
2.53T0.13 3.14T0.01 3.78T0.12 6.3T0.9 
2.70T0.10 3.14T0.02 3.72T0.14 5.8T0.7 

3.13±0.02 3.78T0.08 6.2T0.6 
2.52T0.17 3.16T0.02 3.73T0.15 6.3T0.8 
3.16T0.11 - 6.9T1.0 

3.11T0.03 - 6.4T1.0 

2.51T0.16 2.14T0.03 2.88±0.06 5.2T0.1 
2.25T0.15 2.17T0.02 2.93T0.06 5.2T0.3 
2.36T0.06 2.18T0.01 2.86T0.15 5.1T0.5 

2.21T0.01 2.90T0.06 5.2T0.3 
2.28T0.17 2.18T0.03 2.90T0.11 5.5±0.2 
2.19T0.16 - 6.0T1.1 


“ PBEO results at 380 K from Ref. [38]. 

^ PBEO results at 380 K from Ref. [39]. 

° Empirical potential structure refinement (EPSR) of neutron diffraction data (0.67 M NaCl solutions) of Ref. [12]. 

X-ray scattering experiment of Ref. [10]. 



FIG. 1: The chlorine-oxygen (upper panel) and chlorine- 
hydrogen (lower panel) radial distribution functions, gci-ii{r) 
and 301 - 0 ( 7 '), of the aqueous chloride ion solution obtained 
from theory via the DFT-based AIMD simulations performed 
in this work. 


dated with the first coordination shell were computed by 
integrating each respective g{r) in Fig. 1 up to the first 
minimum via Ua-g = dvrp J dr ga-p{r), wherein p is the 
number density of the respective atom. 


The first well-pronounced feature to note in Fig. 1 is 
the sharp first peak in the gci-oi'i') corresponding to the 
first solvation shell, which is followed by a lower and 
broader second peak representing the second coordina¬ 
tion shell. In general, all of the AIMD trajectories con¬ 
sidered here (generated using the four different XC func¬ 
tionals discussed above) yield consistent RDFs, with each 
predicting that the position of the first maximum is lo¬ 
cated at = 3.14 ± 0.02 A (rg^^n = 2.17 ± 0.02 

A). However, there is a noticeable trend in that PBEO 
and the two vdW-inclusive XC functionals predict val¬ 
ues of and rg'gfjj that are larger by 0.03-0.05 A in 

comparison to PBE. In this regard, all of these values ob¬ 
tained from our simulations are in fairly good agreement 
with the experimental scattering values of 3.11-3.16 A 
for and 2.19 A for (see Table I). Further¬ 
more, we note that our findings for the and rg'gfjj 

quantities at the PBEO level of theory are also in rela¬ 
tively good agreement with the analogous data provided 
by the recent work of Zhang et which studied the 

solvation structure of Cl“ in water using the PBEO XC 
functional (see Table I). With deviations of 0.01 A and 
0.03 A, respectively, these small differences may be at¬ 
tributed to the ~ 50 K temperature difference between 
these two simulations. 

Since the first peaks of the (7ci-o(''’) and 3 ci-h( 7') are 
useful indicators of the relative structuring of the solva¬ 
tion shell surrounding a given Cl“ ion, the dependence 
of these quantities on the underlying XC potential de¬ 
serves further comment. For one, we note that the in¬ 
clusion of vdW/dispersion interactions leads to a reduc¬ 
tion in the structure of the first solvation shell—an ef¬ 
fect that manifests as a simultaneous reduction in the 
intensities coupled with an increase in the positions of 
the hrst maxima in both 301 - 0 ( 7 ') and 3ci-h(7'). At the 
PBE (PBEO) level of theory, the inclusion of non-local 
vdW/dispersion interactions leads to significant reduc¬ 
tions in and hggfjj of 0.14 (0.18) and 0.26 (0.08), 

respectively. Interestingly, the trend observed here in- 













5 


dicates that the vdW/dispersion correction has a sim¬ 
ilar effect on when used in conjunction with ei¬ 

ther a GGA (PBE) or a hybrid (PBEO) XC functional, 
yet the same correction clearly has a more significant 
effect on when used in conjunction with the PBE- 

GGA functional. Since the inclusion of exact exchange 
at the PBEO level of theory weakens the hydrogen bond 
strength with respect to this directly af¬ 

fects the interactions between the Cl~ ion and the neigh¬ 
boring hydrogen atoms in the predominantly hydrogen- 
bonded first solvation shell. Therefore, signifi¬ 

cantly reduced with PBEO compared to PBE (by 0.15), 
whereas the effect of PBEO on marginal 

(by 0.03). In other words, PBEO significantly softens the 
first peak in the ( 7 ci-H(r’), but has little effect on the 
first peak in the gci-o(r)- Hence, the effects of non-local 
vdW/dispersion interactions on significantly 

different when used in conjunction with PBE and PBEO, 
while their effects on quite similar. When com¬ 

pared to PBE, the collective effects of exact exchange and 
vdW/dispersion interactions, i.e., as accomplished with 
the use of the PBEO-l-TS-vdW(SC) XC functional, re¬ 
duce by 0.15 and 0.23, accompanied 

by small increases in and of 0.05 A and 

0.04 A, respectively. From these two observations, it is 
clear that the inclusion of exact exchange and non-local 
vdW/dispersion interactions is effectively weakening the 
interaction between the Cl~ ion and the surrounding wa¬ 
ter molecules in the first solvation shell. 

The effects of vdW/dispersion interactions are also 
clearly evident in the interstitial region between the first 
and second coordination shells in both gc\-o{r) and 
gci-uir). Here, the attractive vdW/dispersion force en¬ 
hances the interaction strength between the Cl~ ion and 
interstitial water molecules and favors more distorted and 
entropically-driven geometrical configurations. Hence, 
we observed a net increase in the population of water 
molecules in the interstitial region in both gci-o(r) and 
Sci-h(''’) with the vdW-inclusive functionals considered 
in this work (see Fig. 1). 

This simultaneous weakening and strengthening of the 
interactions between the Cl~ ion and its surrounding 
water molecules as obtained with vdW-inclusive hybrid 
DFT can be further understood by considering a par¬ 
tial decomposition of the RDFs. In Fig. 2, the contri¬ 
butions to the gci-oif) and (;ci-h(»") arising from the 
4th, 5th, 6th, and 7th neighboring water molecules are 
plotted for each of the XC functionals employed in this 
work. At the PBEO-l-TS-vdW(SC) level of theory, the 
first five oxygen and hydrogen atoms (with only the 4th 
and 5th neighbors shown for clarity) were found to be 
further from the Cl“ ion on average, as compared to the 
corresponding distances obtained with PBE, indicating 
a weakening of the interactions between the CR ion and 
these surrounding water molecules. However, the exact 
opposite trend was found for the 6th and 7th neighbor¬ 
ing oxygen and hydrogen atoms, in which these neighbors 
were found to be much closer to the Cl“ ion in the case 


— PBE - PBE+TS-vdW(SC) - PBEO - PBEO+TS-vdW(SC) 



FIG. 2: Contributions from the 4th (n = 4), 5th (n = 5), 6th 
(n = 6), and 7th (n = 7) nearest neighbors to the Cl-O (left 
panel) and Cl-H (right panel) radial distribution functions, 
pci-o(r) and gci-H(r), of the aqueous chloride ion solution 
obtained from theory via the DFT-based AIMD simulations 
performed in this work. 

of PBEO-l-TS-vdW(SC) when compared to PBE, indicat¬ 
ing a strengthening of the interactions with the Cl“ ion 
and these surrounding water molecules. As a result, the 
collective effects of exact exchange and vdW/dispersion 
interactions tend to weaken the interactions between the 
CR ion and the first five neighbors while simultaneously 
strengthening the interactions with the water molecules 
beyond the 6th neighbors, which results in a larger over¬ 
lap between the first and second coordination shells and 
a net reduction in the structure of the solvation shell sur¬ 
rounding the CR ion. 

Furthermore, the collective effects of exact exchange 
and vdW/dispersion interactions also influence the posi¬ 
tions of the first minima in the gci-oii^) und gci-H(*’), 
quantities which are typically used to indicate the points 
of separation between the first and second coordination 
shells. In this regard, the positions of the first minima in 
these RDFs, i.e., r’c^o ^ci-H; increase by 0.09 and 
0.02 A, respectively, with PBEO-l-TS-vdW(SC) in com¬ 
parison to PBE (see Table I). This finding has an impor¬ 
tant immediate consequence, in that the mean coordina¬ 
tion numbers obtained by integrating the first peak of 
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FIG. 3: Probability distributions of the (a) Cl-0 and (b) Gi¬ 
ll coordination numbers, nci-o and nci-n, calculated by in¬ 
tegrating gci-o(^) and (;ci-h(’') up to their respective first 
minimum. 


these RDFs up to their first respective minima vary con¬ 
siderably among the XC functionals considered herein. 
Here, we found that the mean Cl-0 coordination number, 
nci-Oj increases by 0.8 with both PBE+TS-vdW(SC) 
and PBEO-l-TS-vdW(SC) in comparison to PBE (see Ta¬ 
ble I) and the resultant mean value of nci-o = 6.3 agrees 
reasonably well with recent experimental findings that 
report coordination numbers between 6.4-6.9 in < 2M 
aqueous NaCl solutions.For the mean Cl-H coordi¬ 
nation number, PBE0+TS-vdW(SC) simulations yield a 
value of nci-H = 5-5, which is the largest value among 
all of the XC functionals considered in this work (all of 
which predict a consistent mean value of « 5.2), and is 
in better agreement with the experimental EPSR value 
of 6.0±1.1.^^ We note in passing that these estimates for 
the mean Cl-0 and Cl-H coordination numbers are con¬ 
sistent with previous AIMD simulations performed with 
GCA and hybrid functionals.28-34.36-39 
In order to further elucidate the instantaneous varia¬ 
tions in the geometry of the first coordination shell sur¬ 
rounding the aqueous Cl“ ion, we now analyze the proba¬ 
bility distributions of the nci-o and nci-n coordination 
numbers (see Fig. 3). From this analysis, we found that 
most (~ 40-50%) of snapshots at the PBE level of the¬ 
ory were characterized by 5-fold coordination for both 
nci-O and nci-H (see Fig. 3), indicating that the Cl“ 
ion is predominantly solvated by five water molecules, 
each of which are pointing one hydrogen atom toward 
the Cl“ ion. At the PBE0+TS-vdW(SC) level of theory, 
we instead found that most (?« 35%) of snapshots can be 
characterized by 6-fold Cl-0 coordination and that both 
5- and 6-fold Cl-H coordinations were found to be equally 
most probable (« 35% each). Qualitatively similar obser¬ 
vations can also be made with the PBEO and PBE+TS- 
vdW(SC) XC functionals, in which nci-o is dominated 
by 6-fold coordination, whereas nci-n is dominated by 5- 
fold coordination. This apparent mismatch between the 


most dominant coordination number in the nci-o ^md 
nci-H probability distributions is reflective of the pres¬ 
ence of more disordered hydrogen-bond structures in the 
first solvation shell. 

To graphically depict the differential coordination 
numbers found in the aqueous Cl“ ion solution, we have 
computed the spatial density functions (SDF) within the 
first solvation shell of the Cl“ ion based on structures 
(configurations) obtained at the PBE and PBEO+TS- 
vdW(SC) levels of theory. In Fig. 4, we have plotted 
SDFs for configurations in which the Cl“ ion was coor¬ 
dinated with 4, 5, 6, and 7 water molecules, respectively, 
averaged over all structures characterized by these coor¬ 
dination numbers. For the case of 4-fold coordination, 
we found that distorted tetrahedral structures are domi¬ 
nant over square planar structures in both XC function¬ 
als (see Fig. 4(a) and Fig. 4(e)). In both XC function¬ 
als, the 5-fold coordinated aqueous Cl“ ion complexes 
are predominantly comprised of distorted square pyra¬ 
midal structures instead of trigonal bipyramidal struc¬ 
tures; however, the SDF of the oxygen and hydrogen 
atoms obtained with PBEO-l-TS-vdW(SC) (see Fig. 4(f)) 
were found to be more delocalized than that obtained 
with PBE (see Fig. 4(b)). For the case of 6-fold coor¬ 
dination, the predominant form found with PBE resem¬ 
bles an octahedral arrangement (see Fig. 4(c)), whereas 
the trigonal prism structure was more representative of 
the PBEO-l-TS-vdW(SC) configurations (see Fig. 4(g)). 
The 7-fold coordination case again significantly differs 
among these distinct XC functionals; here the predomi¬ 
nant shape found at the PBE level of theory resembles a 
capped octahedron (see Fig. 4(d)), while the more com¬ 
mon pentagonal bipyramidal arrangement dominates at 
the PBE0+TS-vdW(SC) level of theory (see Fig. 4(h)). 
In all of the SDF plots in Fig.. 4, the oxygen and hy¬ 
drogen atoms were found to be more delocalized on av¬ 
erage in the PBEO-l-TS-vdW(SC) configurations. This 
observation, in conjunction with the earlier finding that 
a majority of the configurations posses 6- and 7-fold Cl- 
O coordination, again exemplifies the fact that the first 
solvation shell surrounding the Cl“ ion is significantly 
more disordered at the vdW-inclusive hybrid PBEO-t-TS- 
vdW(SC) level of theory. 

In addition to the detailed coordination structures con¬ 
sidered above, the aforementioned disorder in the first 
coordination shell surrounding the Cl“ ion can also be 
characterized by analyzing the O-Cl-0 angular distribu¬ 
tion function. In this regard, O-Cl-0 angular distribution 
functions were computed for each of the XC functionals 
employed herein by considering only the oxygen atoms 
which reside within the first coordination shell surround¬ 
ing the Cl“ ion (defined by the radial cutoff distance 
given by As seen in Fig. 5, the PBE XC func¬ 

tional predicts a maximum at 0 « 90°, a finding which 
is consistent with the fact that a majority of the con¬ 
figurations at the PBE level are characterized by 5-fold 
Cl-0 coordination (see Fig. 3) and predominantly square 
pyramidal structures (see Fig. 4(b)). On the other hand. 
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FIG. 4: Spatial density functions (SDF) of the first solvation shell surrounding the aqueous Cl“ ion characterized by 4-, 5-, 
6-, and 7-fold coordination numbers, respectively (shown from left to right). The top ((a)-(d)) and bottom ((e)-(h)) panels 
correspond to structures (conhgurations) obtained from PBE and PBEO-l-TS-vdW(SC) AIMD simulations, respectively. The 
Cl“ ion, oxygen atoms, and hydrogen atoms are depicted by green, red, and cyan spheres, respectively. 


PBEO and both vdW-inclusive XC functionals predict 
maxima at 0 ss 75°, which is consistent with the fact 
that a majority of the configurations at these levels of 
theory are characterized by 6-fold Cl-0 coordination (see 
Fig. 3) and predominantly distorted trigonal prism struc¬ 
tures (see Fig. 4(g)). Hence, this analysis of the O-Cl-0 
angular distribution functions is again strongly indicative 
of the net increase in the amount of structural disorder 
present in the first solvation shell surrounding the Cl~ ion 
when both exact exchange and non-local vdW/dispersion 
interactions are accounted for in the underlying XC po¬ 
tential. 


B. Electronic Structure of the Aqueous Cl Ion 

Beyond the solvation structures of aqueous ionic so¬ 
lutions, which can be directly compared to diffraction 
experiments, the electronic properties of aqueous ionic 
solutions, such as the band structure of the solution and 
the energy levels of the solvated ion, are also important 
quantities that can be probed by photoemission (PES) 
and X-ray absorption spectroscopy (XAS). In this regard, 
a recent PES experiments^ has measured electron bind¬ 
ing energies in aqueous chloride ion solutions, in which 
the position of the highest occupied molecular orbital 
(HOMO) of the solvated chloride ion was found to be 
5E = 1.25-1.50 eV above the valence band maximum 
(VBM) of liquid water. In order to theoretically study 
the electronic structure properties of the aqueous Cl“ 



FIG. 5: Angular distribution functions, O-Gl-O, for each of 
the XC functionals employed herein computed by considering 
only the oxygen atoms which reside within the first coordina¬ 
tion shell surrounding the Gl“ ion (as defined by r^l^Q, the 
position of the first minimum in the (;ci-o(r)). 


ion solution, we have performed DFT-based ground state 
wavefunction optimizations based on equilibrated struc¬ 
tures (configurations) from various AIMD trajectories. 
In particular, we have also adopted several different XC 
approximations, ranging from semi-local GCAs (PBE) 
to hybrid functionals (PBEO and BHLYP) in these band 
structure calculations. 

We begin this analysis by first considering the band 
gaps (HOMO-LUMO gaps) of liquid water and the aque- 












TABLE II: Electronic structure (ES) properties computed 
with the PBE, PBEO, and BHLYP XC functionals based on 
structures (configurations) generated with AIMD simulations 
at the PBE, PBE+TS-vdW(SC), PBEO, and PBEO+TS- 
vdW(SC) levels of theory. and are the computed 

band gaps (in eV) for liquid water and the aqueous Cl“ ion 
solution, respectively. dE is the energetic difference between 
the HOMO of the Cl“ ion and the valence band maximum 
(VBM) of liquid water. 


ES Level//AIMD Level 

pti2^ 

Ef 

5E 

PBE//PBE 

- 

4.22 

-0.24 

PBE0//PBE 

- 

6.68 

0.13 

BHLYP//PBE 

- 

8.81 

0.72 

PBE//PBE+TS-vdW(SC) 

- 

4.16 

-0.11 

PBEO //PBE-ETS-vdW(SC) 

- 

6.53 

0.27 

BHLYP//PBE+TS-vdW(SC) 

- 

8.61 

0.88 

PBE//PBE0 

- 

4.19 

-0.06 

PBE0//PBE0 

- 

6.55 

0.35 

BHLYP//PBE0 

- 

8.66 

0.95 

PBE //PBEO-ETS-vdW(SC) 

4.42 

4.20 

-0.06 

PBE0//PBE0+TS-vdW(SC) 

7.04 

6.57 

0.34 

BHLYP//PBEO-ETS-vdW(SC) 

9.70 

8.62 

0.92 


ous Cl“ ion (see Table II). Here, it is well-known that the 
accuracy of band structure calculations depends crucially 
on the XC approximation employed in the quantum me¬ 
chanical description of the electrons. In this regard, the 
inherent delocalization error and the resulting incorrect 
convex behavior characteristic of semi-local GGA func¬ 
tional approximations leads to large underestimates of 
band gaps and incorrect predictions for other features in 
the band structure.Gonsistent with previous DFT cal¬ 
culations of liquid water employing the PBE functional,^® 
the estimated band gap was found to be ~ 4.4 eV in this 
work, which is underestimated by 4.3 eV from the best 
experimental estimates.Incorporation of a fraction of 
exact exchange {E^) is an established approach for rem¬ 
edying this problem to some extent,yielding band 
gaps of approximately 7.0 and 9.7 eV for liquid water us¬ 
ing the PBEO (25% E^x) and BHLYP (50% Exx) hybrid 
XC functionals, respectively (see Table H). Hence, the in¬ 
clusion of Exx reduces the delocalization (self-interaction) 
error and leads to a marked improvement in the predicted 
band gaps of liquid water. In this regard, it should also 
be noted that the computed band gaps are relatively in¬ 
sensitive to the molecular structures (configurations) ob¬ 
tained from AIMD trajectories generated by various XC 
functional approximations.®® 

We now focus our attention on the electronic proper¬ 
ties of the aqueous Cl~ ion solution. In particular, 5E^ 
which is defined as the relative energetic difference be¬ 
tween the solvated Cl“ ion and the VBM of liquid water 
will play a key role in our analysis. In this regard, it 
was previously found that the energy level of the sol¬ 
vated Cl“ ion is quite sensitive to both the underlying 
molecular structure as well as the choice of the electronic 


structure method employed during ground state wave- 
function optimization.®® Therefore, a systematic study 
that differentiates these two effects is necessary and will 
be addressed in this work. 



Energy Eigenvalues (eV) 


FIG. 6: Total density of states (DOS in red) and real-space 
projected density of states (PDOS in blue) for the aqueous 
Cl“ ion solution computed using the hybrid PBEO XC func¬ 
tional based on structures (configurations) generated from 
AIMD trajectories at the (a) PBE, (b) PBE-tTS-vdW(SC), 
(c) PBEO, and (d) PBEO-ETS-vdW(SC) levels of theory. The 
DOS and PDOS above were obtained from electronic struc¬ 
ture calculations on 200 configurations taken from each AIMD 
trajectory. 


To accomplish this goal, we first performed band 
structure calculations with the hybrid PBEO XG func¬ 
tional based on structures (configurations) generated 
with AIMD simulations at the PBE, PBE-l-TS-vdW(SC), 
PBEO, and PBEO-l-TS-vdW(SC) levels of theory. With 
a self-consistent solution to the non-linear Kohn-Sham 
(KS) equations, we projected the optimized ground-state 
KS eigenfunctions onto the solvated Gl“ ion with a ra¬ 
dius of i? = 1.5 The resulting total density of states 
(DOS) and real-space projected DOS (PDOS) computed 
at the PBEO level of theory are plotted in Fig. 6 , from 
which it is evident that the total DOS of the aqueous Gl“ 
ion solution is dominated by the features belonging to the 
valence electrons of liquid water. In this regard, the oxy¬ 
gen 2p states are located well above (> 10 eV) the oxy¬ 
gen semi-core 2 s states, and can be further decomposed 
into three features, corresponding to the I 62 , 3ai, and 
Ibi bands representing the covalently bonded and non- 
bonded (lone pair) electrons. As far as the PDOS of the 
solvated Gl“ ion is concerned, the 3s states of the Gl“ ion 
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were found to be between the 2s and 2p bands of liquid 
water while the 3p states of the Cl“ ion were found to be 
located at the edge of the VBM of liquid water—a picture 
which is consistent with the PES measurements in aque¬ 
ous Cl“ ionic solutions.Moreover, it can be seen that 
the total DOS computed based on molecular configura¬ 
tions generated from different AIMD trajectories remain 
very similar in terms of the peak positions and spectral 
distributions. This finding is indicative of the fact that 
the band structure of liquid water is rather insensitive to 
the underlying molecular configurations (ie., to whether 
or not Exx and/or non-local vdW/dispersion interactions 
are accounted for in the underlying XC potential). 

In sharp contrast, the relative position of the orbital 
levels associated with the solvated Cl“ ion with respect 
to the VBM of liquid water shows a much stronger de¬ 
pendence on the underlying molecular configurations. As 
plotted in Fig. 6, the PDOS computed at the PBEO level 
of theory based on molecular configurations generated 
from the PBE-AIMD trajectory is mostly located below 
the VBM of liquid water, yielding a value of 6E = 0.13 
eV (see Table II). With exact exchange and/or non-local 
vdW/dispersion interactions accounted for in the under¬ 
lying potential to generate the AIMD trajectories, the 
distribution of the PDOS becomes more prominantly sep¬ 
arated from the VBM of liquid water. As a result, the 
computed SE increases in the direction of the experimen¬ 
tal value of SE = 1.25-1.50 eV as shown in Table II, a 
result that we attribute to the weakened interaction be¬ 
tween the Cl“ ion and its surrounding water molecules 
with the use of vdW-inclusive hybrid XC functionals. As 
discussed above in Sec. Ill A, the collective effects of E^x 
and non-local vdW/dispersion interactions lead to a net 
decrease in the strength of the hydrogen bonds exist¬ 
ing between the lone pair electrons of the Cl“ ion and 
the water molecules residing in the first solvation shell. 
With the positions of the first maximum in gci- oir) and 
Sci-h(»'), he., and slightly increased, hy¬ 

bridization between the Cl“ 3p and the oxygen 2p or¬ 
bitals (located at the edge of the liquid water valence 
band) becomes less favorable. We note in passing that 
this explanation is consistent with the increased SE de¬ 
picted in Fig. 6. 

From the above discussion, one can see that an ac¬ 
curate prediction of the electronic properties of the sol¬ 
vated Cl“ ion crucially depends on the description of 
the interaction between the p electrons of Cl and O. In 
this regard, the more accurate prediction of the hydrogen 
bond tends to discourage hybridization among these p 
orbitals, thereby increasing SE towards the experimental 
value. This trend is even more pronounced when different 
XC approximations are utilized in the theoretical treat¬ 
ment of the electronic structure. As shown in Fig. 7, we 
now present the DOS and PDOS of the aqueous CD ion 
solution based on ground-state wavefunction optimiza¬ 
tions performed at the PBE, PBEO, and BHLYP levels 
of theory. In this case, the underlying molecular struc¬ 
tures or configurations are kept fixed and were taken from 



Energy Eigenvalues (eV) 


FIG. 7: Total density of states (DOS in red) and real-space 
projected density of states (PDOS in blue) for the aqueous 
CD ion solution computed using the (a) PBE, (b) PBEO, and 
(c) BHLYP XC functionals based on structures (configura¬ 
tions) generated from an AIMD trajectory at the PBEO+TS- 
vdW(SC) level of theory. The DOS and PDOS above were 
obtained from electronic structure calculations on 200 config¬ 
urations taken from the AIMD trajectory. 

an AIMD trajectory employing the vdW-inclusive hybrid 
PBEO-l-TS-vdW(SC) XC potential. This choice is moti¬ 
vated by our recent study in which we found that the 
PBEO-l-TS-vdW(SC) XC potential yielded a very accu¬ 
rate prediction of the microscopic structure of ambient 
liquid water that is in quantitative agreement with re¬ 
spect to the available scattering experiment data.^°^ 

As clearly seen from Fig. 7, the PDOS of the sol¬ 
vated CD ion is significantly affected by the choice of 
DFT functional employed in the electronic structure cal¬ 
culation. At the semi-local PBE-GGA level of theory, 
the computed value of SE = —0.06 eV is negative, 
which is opposite in sign compared to the experimen¬ 
tal value; however, as the fraction of Exx is increased 
from 25% (PBEO) to 50% (BHLYP), the SE steadily 
increases to 0.34 eV and 0.92 eV, respectively, and is 
rapidly approaching the experimental value. Again, we 
attribute this trend to the removal of the deleterious self¬ 
interaction error that is present in semi-local XG func¬ 
tionals via the admixture of exact exchange in the hybrid 
PBEO and BHLYP XG functionals. As such, the elec¬ 
tronic structure predicted by semi-local GGA function¬ 
als (such as PBE) artificially favors the interaction of the 
GD ion and its surrounding water molecules by facilitat¬ 
ing a greater extent of p orbital hybridization. This effect 
is greatly reduced at the hybrid DFT level of theory— 
at the BHLYP level, the eigenfunctions of the GD ion 
HOMO were found to have a non-negligible amplitude 
located on water molecules in the first coordination shell 
in only « 5% of the configurations. However, we note in 
passing that this fraction increases to over 70% when the 
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semi-local PBE-GGA XC functional is used instead (see 
Fig. 8). 



FIG. 8: The triply degenerate gas-phase (left) and condensed- 
phase (right) 3p orbitals for the Cl“ ion. 

These findings reported above are consistent with the 
previous study by Zhang et al.^^ in which it was found 
that the use of hybrid XC functionals in generating 
both the electronic and molecular structure was essen¬ 
tial for a qualitatively correct prediction of 5E. The 
current study further indicates that accounting for non¬ 
local vdW/dispersion interactions in the generation of 
the underlying molecular structure is also important for 
obtaining a quantitatively correct theoretical prediction 
of this fundamental quantity in aqueous ionic solutions. 
By admixing 25% of Fxx, the hybrid PBEO XC func¬ 
tional has successfully generated an electronic structure 
that is qualitatively consistent with the available PES ex¬ 
periments. However, 25% of might not be sufficient 
to completely eliminate the effects of self-interaction er¬ 
ror on which can serve as an explanation for the 

remaining discrepancy that exists between theory and 
experiment. In addition, it should also be noted that an 
explicit treatment of nuclear quantum effects (NQE), as 
provided by the Feynman discretized path integral tech¬ 
nique, is expected to further decrease the strength of the 
interaction between the Cl“ ion and water—as a result, 
a further increase in 6E towards the experimental value 
is anticipated with such a theoretical treatment of this 
fundamental aqueous ionic solution. 


IV. CONCLUSIONS 

In this work, we have systematically studied the sol¬ 
vation and electronic structure of the aqueous chloride 
ion solution using DFT-based AIMD which incorporates 
exact exchange (Exx) and non-local vdW/dispersion in¬ 
teractions in the underlying XC potential. From an anal¬ 
ysis of gci-oif) and gci-H (''’), we find that the inclusion 
of exact exchange and non-local vdW/dispersion interac¬ 
tions via the use of the PBEO-l-TS-vdW(SC) XC func¬ 
tional effectively weakens the interactions between the 
Cl“ ion and the surrounding water molecules residing in 
the first solvation shell and yields a mean Cl-0 coordina¬ 
tion number, nci-O = 6.3±0.8, that is in excellent agree¬ 
ment with the available experimental data. Moreover, 


we found that most configurations at the vdW-inclusive 
hybrid PBEO-l-TS-vdW(SC) level of theory can be char¬ 
acterized by 6-fold Cl-0 coordination and are predomi¬ 
nantly comprised of distorted trigonal prism structures— 
findings which are strongly indicative of a significantly 
disordered first solvation shell surrounding the Cl“ ion. 

By performing a systematic series of band structure 
calculations on molecular structures (configurations) gen¬ 
erated from AIMD simulations with various differing XC 
potentials, we were able to selectively isolate the effects 
arising from the underlying molecular structure as well 
as the choice of the electronic structure method on the 
energy levels of the solvated Cl“ ion relative to the VBM 
of liquid water {5E). In doing so, we found that while 
the band structure of liquid water is rather insensitive 
to the underlying molecular configurations, the relative 
positions of the orbital levels associated with the sol¬ 
vated Cl“ ion with respect to the VBM of liquid wa¬ 
ter strongly depend on the underlying molecular config¬ 
urations. Here, the effective weakening of the hydrogen 
bonds existing between the lone pair electrons of the Cl“ 
ion and the water molecules residing in the first solva¬ 
tion shell due to the collective effects of Fxx and non¬ 
local vdW/dispersion interactions leads to a decrease in 
the degree of hybridization between the Cl“ 3p and the 
oxygen 2p orbitals located at the edge of the liquid water 
valence band. In the same breath, the PDOS of the sol¬ 
vated Cl“ ion is also significantly affected by the choice of 
DFT functional employed in the electronic structure cal¬ 
culation; as the fraction of E^x was increased from 25% 
(PBEO) to 50% (BHLYP), SE steadily increased from 
0.34 eV to 0.92 eV, thereby rapidly approaching the ex¬ 
perimental value of 1.25-1.50 eV. 

As a hnal set of remarks, we note that the quantitative 
description of the microscopic structure of the aqueous 
Cl“ ion solution may change by further refinement of the 
underlying XC functional approximation, which may be 
accomplished by reducing the self-interaction error via 
fine-tuning the exchange correction,and/or includ¬ 
ing a better description of vdW/dispersion interactions 
via the inclusion of beyond-pairwise interactions, as pro¬ 
vided, for example, by the recently proposed many-body 
dispersion (MBD) scheme.As mentioned above, 
structural refinements are also expected to result from 
a proper treatment of nuclear quantum effects via the 
Feynman discretized path-integral approach, which is ex¬ 
pected to further weaken the interactions between the 
solvated Cl“ ion and its aqueous surroundings; the use 
of such configurations in conjunction with a higher-level 
electronic structure theory, such as Hedin’s GW self¬ 
energy approximation^^^ and/or the random phase ap¬ 
proximation (RPA),^^^“^^® should lead to even better 
agreement with PES experiments on the relative ener¬ 
getics among the solvated Cl“ orbitals and the VBM of 
liquid water. 
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